The Theory and Simulation of the 21 -cm 
Background from the Epoch of Reionization 

Paul R. Shapiro*, Ilian T. Iliev^ Garrelt Mellema**, Ue-Li Pen* and Hugh 

Merz § 

* Department of Astronomy, University of Texas, Austin, TX 78712-1083, U.S.A. 
^Institute of Theoretical Physics, University Zurich, Winterthurerstrasse 190, CH-8057 Zurich, 

Switzerland 

"Stockholm Observatory, AlbaNova University Center, Stockholm University, SE-106 91 

Stockholm, Sweden 

^Canadian Institute for Theoretical Astrophysics, University of Toronto, 60 St. George Street, 

Toronto, ON M5S 3H8, Canada 
^Department of Physics & Astronomy, University of Waterloo, 200 University Avenue West, 

Waterloo, ON N2L 3G1, Canada 

Abstract. The redshifted 21 -cm line of distant neutral H atoms provides a probe of the cosmic 
"dark ages" and the epoch of reionization ("EOR") which ended them, within the first billion years 
of cosmic time. The radio continuum produced by this redshifted line can be seen in absorption or 
emission against the cosmic microwave background ("CMB") at meterwaves, yielding information 
about the thermal and ionization history of the universe and the primordial density perturbation 
spectrum that led to galaxy and large-scale structure formation. Observing this 21 -cm background 
is a great challenge, as it is necessary to detect a diffuse signal at a brightness temperature that differs 
from that of the CMB at millikelvin levels and distinguish this from foreground continuum sources. 
A new generation of low-frequency radio arrays is currently under development to search for this 
background. Accurate theoretical predictions of the spectrum and anisotropy of this background, 
necessary to guide and interpret future observations, are also quite challenging. Toward this end, it 
is necessary to model the inhomogeneous reionization of the intergalactic medium and determine 
the spin temperature of the 21 -cm transition and its variations in time and space as it decouples from 
the temperature of the CMB. In my talk, I summarized some of the theoretical progress in this area. 
Here, I will focus on just a few of the predictions for the 21-cm background from the EOR, based 
on our newest, large-scale simulations of patchy reionization. These simulations are the first with 
enough N-body particles (from 5 to 29 billion) and radiative transfer rays to resolve the formation 
of and trace the ionizing radiation from each of the millions of dwarf galaxies believed responsible 
for reionization, down to 10 8 M Q , in a cubic volume large enough (90 and 163 comoving Mpc on a 
side) to make meaningful statistical predictions of the fluctuating 21-cm background. 
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INTRODUCTION 

The redshifted 21-cm line of neutral hydrogen atoms can probe the thermal and ioniza- 
tion history of the universe from the cosmic "dark ages" through the EOR and the density 
perturbations responsible for structure formation (see, e.g., yj] for review and refs). It 
is generally believed that, during the EOR, the spin temperature, 7^ , of this hyperfine 
transition was decoupled from the CMB temperature, 7cmb> by the Wouthuy sen-Field 



mechanism, whereby UV photons at the Ly a transition frequency are resonantly scat- 
tered, and the excited atoms then decay to a different hyperfine level of the ground state. 
This radiative pumping drives T s toward the gas kinetic temperature, Tk- The 21-cm 
line is seen in emission or absorption against the CMB if 7^ is greater than or less than 
7cmb> respectively. If, as is also generally believed, other radiation is also present (e.g. 
X-rays) to heat the neutral IGM so that Tk 3> 7cmb> men tne 21-cm differential bright- 
ness temperature depends only upon the neutral gas density and peculiar velocity fields, 
independent of the value of Tk or T s . Henceforth, we shall limit our discussion of the 
21-cm background from the EOR to this regime. 

Reionization is generally believed to be the outcome of the release of ionizing radia- 
tion by galaxies undergoing star formation (see, e.g., [0, Hi for recent reviews). Current 
theory suggests that the galaxies responsible for most of this radiation are dwarf galaxies 
more massive than about 1O 8 M . While dark matter dominates the gravitational forces 
which cause this structure formation, ordinary atomic matter must join the dark mat- 
ter in making galaxies for star formation to be possible. Once the atomic gas in the 
intergalactic medium ("IGM") in some region is heated by reionization, however, gas 
pressure opposes gravitational collapse, and, thereafter, the smallest galaxies form with- 
out atomic matter and cannot make stars. The minimum mass of star- forming galaxies 
in such regions is about 10 9 solar masses. This suppression of star-forming galaxies of 
lower mass inside ionized regions is sometimes referred to as "Jeans-mass filtering" [0]. 

Due to its complexity, reionization is best studied through large-scale numerical 
simulations D5ilaLZil8|], combining the numerical challenges of Af-body simulation and 
radiative transfer. The tiny galaxies which are the dominant contributors of ionizing 
radiation must be resolved in very large cosmological volumes, large enough to contain 
billions of times more total mass than one dwarf galaxy and up to tens of millions of 
such galaxies, in order to derive their numbers and clustering properties correctly. The 
expansion of the intergalactic ionization fronts created by all these millions of galaxies 
must, in turn, be tracked with a 3D radiative transfer method which also solves non- 
equilibrium ionization rate equations. The combination of all these requirements makes 
this problem formidable. 

Our original simulations [5s] resolved the formation of all galaxies more massive than 
about one billion solar masses, those which are expected to form stars even after reion- 
ization has heated their environment. These simulations of a comoving, cosmologically- 
expanding cubic volume (150 Mpc) 3 enabled us to make the first statistically meaning- 
ful predictions of the impact of reionization on observables like the 21-cm background 
10], the CMB anisotropy fTToh . and the observability of early galaxies as Lyman a line- 
emission sources yjj], essential predictions to guide and interpret current and future 
observational programs |fl3| . With the next generation of simulations described here, 
we shall explore the role of less massive dwarf galaxies. These are important sources 
of ionizing radiation if they form before their neighborhood is reionized but are pre- 
vented from being sources if they form after it is reionized. Our previous results showed 
that the inclusion of these sources and of their suppression changes the outcome of 
reionization substantially [13]. Simulations of volume (50 Mpc) 3 showed that reioniza- 
tion was "self-regulated," i.e. the more abundant low-mass halos initially dominated but 
later choked themselves off as their formation was suppressed in ionized regions. The 
next step, described here, is to push these simulations of "self-regulated" reionization to 



larger volumes, large-enough to make meaningful statistical predictions of observables 
like the 21 -cm background. 



NEW, LARGE-SCALE SIMULATIONS OF SELF-REGULATED 

REIONIZATION 

We have developed a new, cosmological Af-body code more powerful than the one we 
used for our previous simulations jl4 \. CubeP 3 M, a massively -parallel P 3 M (particle- 
particle-particle-mesh) lfl5ll code, is a successor to the previous code, PMFAST lfl6ll . 
Its force resolution is significantly better than that of PMFAST, by inclusion of the 
short-range correction from direct particle-particle (PP) forces. With this improved force 
resolution, dark mater halos can be identified with as few as 20 particles per halo, to 
locate all the galaxies which are possible sources of ionizing starlight, with masses above 
this minimum. 

These Af-body results provide the evolving density field of the IGM and the loca- 
tion and mass of all the halo sources, as input to a separate radiative transfer simula- 
tion of inhomogeneous reionization. The latter simulation is performed by our C 2 Ray 
(Conservative, Causal Ray-Tracing) code, a grid-based, ray-tracing, radiative transfer 



and nonequilibrium chemistry code, described in | |17|1 . The ionizing radiation is ray- 
traced from every source to every grid cell using a method of short characteristics. The 
code is explicitly photon-conserving in both space and time, which ensures an accurate 
tracking of ionization fronts, independent of the spatial and time resolution, even for 
grid cells which are optically thick to ionizing photons and time steps long compared to 
the ionization time of the atoms, with corresponding great gains in efficiency. The code 
has been tested against analytical solutions QJ/ZU and, in direct comparison with other ra- 
diative transfer methods, on a standardized set of benchmark problems II 1 811 . C 2 Ray has 
also been used in the first, 3D radiation-hydrodynamical simulations of an H II region 



in a turbulent molecular cloud [19]. While our basic methodology remains essentially as 
described in II 1711 . our C 2 Ray code is now parallelized for distributed-memory machines 

Q. 

With CubeP 3 M, we simulated the ACDM universe with 1728 3 and 3072 3 Af-body par- 
ticles of mass 5 x 1O 6 M , in comoving box sizes of 91 Mpc and 163 Mpc, respectively. 
At 29 billion particles, the latter is the largest N-body simulation of the formation of 
early cosmological structure to date. This mass resolution allows us to resolve all ha- 
los above 10 8 M Q , roughly the minimum mass of halos which can radiatively cool by 
hydrogen atomic-line excitation and efficiently form stars. These halos are separated 
into halos above the Jeans mass of the ionized IGM (~ 10 9 M Q ), which are not affected 
by reionization, and lower-mass halos, whose star formation is suppressed if the halo 
formed inside an intergalactic H II region. 

Figured] (a) shows the halo mass functions at z = 6 for three of our N-body simula- 
tions. These Af-body simulations, i.e. the halo catalogues and density field, are the basis 
for our C 2 Ray radiative transfer simulations of reionization. All halos are assumed to 
host galaxies, and thus, ionizing sources, unless they form inside an HII region with a 
mass below the threshold for "Jeans-mass filtering". To each galaxy, we assign an ioniz- 
ing photon luminosity proportional to the halo mass, with an efficiency, fy, which gives 
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FIGURE 1. (a)(left) Halo mass functions at z = 6 from A^-body simulations of ACDM with box 
sizes 91 Mpc (blue triangles), 103 Mpc (green squares) and 163 Mpc (red circles). The three mass 
functions agree with each other, except for the largest halos, which are subject to cosmic variance 
for the smaller volume simulations. The lines show the predicted mass functions based on the well- 
known analytical models of Press-Schechter ([20]; dotted) and Sheth-Tormen ([21]; dashed), (b) (right) 
(top panels) Total number of halos (dashed) and number of active sources (i.e. total number minus the 
suppressed ones; solid); and (bottom panels) collapsed mass fraction in high-mass (red, solid) and low- 
mass (i.e. suppressible) sources (blue, dashed) (insets show the same in log scale) vs. cosmic redshift for 
simulations fl00_250S_91Mpc (left) and f50_250S_163Mpc (right). 



the number of photons emitted per atom over the lifetime of the sources (here assumed to 
be 17.6 and 1 1.5 Myr for our 91 Mpc and 163 Mpc simulations, respectively). The sim- 
ulations are denoted fl00_250S_91Mpc_432 and f50_250S_163Mpc_384, where the 
first two numbers are the values of fy for the high- and low-mass sources, respectively, 
"S" indicates that source suppression due to Jeans-mass filtering was included, the third 
is the box size and the last is the radiative transfer grid resolution per dimension. 

Figure |TJb) shows the evolution of the total source numbers, active source numbers 
and the fractions of the mean mass density collapsed into the two types of sources, 
those too massive to be suppressed by reionization and those which are sources only 
if they form in a neutral patch of the IGM. Figure [2] (and accompanying animationj) 
shows images of the gas density field and (active) source halos overlayed with the 
ionized fraction field derived from our radiative transfer simulations. These indicate the 
complex, evolving geometry and topology of the network of ionized bubbles and neutral 
patches, which we must calculate in detail in order to make predictions of the fluctuating 
21 -cm background from this EOR. 

The evolution of the mean ionized fractions of the IGM, mass-weighted (x m ) and 
volume weighted (x v ), is shown in Figure [31 including a comparison with those from 
our previous simulations of "self-regulated" reionization in a smaller box (50 Mpc), de- 
scribed in [13], and those from our previous large-box (143 Mpc) simulations with halo 
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FIGURE 2. Spatial slices of the ionized and neutral gas density from our radiative transfer simulations 
(a)(left) boxsize 163 Mpc at redshift z = 11.6 and (b)(right) boxsize 91 Mpc at z = 11.9, both at 
box-averaged ionized fraction by mass x m — 0.30. Shown are the density field (green) overlayed with 
the ionized fraction (red/orange/yellow) and the cells containing sources (dark/blue). There is a good 
agreement between the two simulations in the typical sizes of the locally-percolated ionization bubbles. 
[A "Quicktime" animation of the simulation on the right is included here.] 



masses above 1O 9 M only (i.e. no suppressible, low-mass halos) Ill2ll . In all cases, reion- 
ization proceeds "inside-out" - i.e. halos form in high-density regions, clustered around 
density peaks and ionize these environments first, so x m > x v . For the new, self-regulated 
case shown there, x m = 0. 1, 0.5, and 0.99 at z = 13.9, 10.7, and 7.9, respectively. As de- 
scribed in |fl3|l . the effect of including the lower-mass (suppressible) halos is to enable 
reionization to start earlier, since the low-mass halos are more abundant earlier. The end 
of reionization is not significantly affected, however, since it is dominated by the expo- 
nential rise of the abundance of the higher-mass halos, which are unaffected by Jeans- 
mass filtering of the IGM, while the abundance of low-mass halo sources saturates as 
the ionized volume grows. 

Predictions of the fluctuating 21 -cm brightness temperature, dTy, from our 91 Mpc 
box (Adb ox ~ 35' at z ~ 8) simulations are shown in Figure HI which maps the sky along 
the line of sight in position - redshift (i.e. angle-frequency) space. The rms fluctuations 
shown in Figures [5ta) and (b) for this case and the 163 box simulations, too, vary 
with frequency as the ionization and density fluctuations evolve in redshift. There is a 
broad but distinct peak in the frequency range Av ~ 120 — 140 Hz (i.e. x m m 0.5 — 0.8). 
Figure |3c) also shows the 2-D angular power spectrum of 57^, when x m = 0.3, which 
peaks at/^3 — 4xl0 4 . Qualitatively similar features have been seen before from 
EOR simulations, but these are the first simulations which are not only of volumes 
large enough to make statistically reliable predictions, but also of mass-resolution high 
enough to take explicit account of the self-regulated nature of reionization. 




FIGURE 3. Mean ionized fractions, x m (mass-weighted; dashed) and x v (volume-weighted; solid) vs. 
redshift for three simulations: new simulations of box size 91 Mpc (fl00_250S_64Mpc_WMAP 3+; 
red; N.B. size label "64" here means "64/h", where h = 0.7); previous "self-regulated reionization" 
simulations of box size 50 Mpc (f250_250S_35Mpc_WMAP3; magenta) from 111311 : previous large- 
volume simulations without low-mass halos, of box size 143 Mpc (f250_100Mpc_WMAP3; blue) from 

m. 
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FIGURE 5. 21 -cm brightness temperature fluctuations (a) (left) rms fluctuations {8T 2 ) l l 2 [mK] ver- 
sus mean ionized fraction x m , for 91 Mpc (solid) and 163 Mpc (dotted) simulations, if observed with a 
beamwidth A0b eam = 3' and bandwidth Av Q b s = 0.2 MHz; (b) (middle) mean differential brightness tem- 
perature 8 T[, [mK] (top) and rms (8T 2 ) 1 / 2 [mK] (middle) versus observed frequency, and, for comparison, 
8Tb{n 2 } x l 2 / (n), to show that signal is enhanced relative to matter fluctuations, by reionization patchiness, 
and evolves differently; (c) (right) 2-D angular power spectrum of 8Tb vs. /, for epoch at which x m = 0.3. 
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